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A substitute for the singular Green kernel 
in the Newtonian potential of celestial bodies 
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ABSTRACT 



The "point mass singularity" inherent in Newton's law for gravitation represents a major difficulty in accurately 
determining the potential and forces inside continuous bodies. Here we report a simple and efficient analytical method 
to bypass the singular Green kernel l/|r — r'\ inside the source without altering the nature of the interaction. We 
build an equivalent kernel made up of a "cool kernel", which is fully regular (and contains the long-range —GM/r 
asymptotic behavior), and the gradient of a "hyperkernel", which is also regular. Compared to the initial kernel, these 
two components are easily integrated over the source volume using standard numerical techniques. The demonstration 
is presented for three-dimensional distributions in cylindrical coordinates, which are well-suited to describing rotating 
bodies (stars, discs, asteroids, etc.) as commonly found in the Universe. An example of implementation is given. The case 
of axial symmetry is treated in detail, and the accuracy is checked by considering an exact potential/surface density 
pair corresponding to a flat circular disc. This framework provides new tools to keep or even improve the physical 
realism of models and simulations of self-gravitating systems, and represents, for some of them, a conclusive alternative 
to softened gravity. 
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1. Introduction 

As a direct consequence of Newton's law for gravitation 
()Newtonlll760t lKellogd[T929[) . the potential of any contin- 
uous distribution of matter inside a volume V at a point 
P(r) of space is given by 



V'(r) 



-G 



p{r')dh 



(1) 



where p{r') is the mass density at P'(r') g V, and (Pt is the 
elementary volunicQ. In general, this is a three-dimensional 
(3D) converging integral. The presence of the Green ker- 
nel l/\r — r'\ is known to represent a difficulty in cal- 
culating -0 everywhere inside and very close to V since 
this function diverges as r — >■ r'. This singularity is clas- 
sically avoided by converting the Green function into a. n 
infinite series (e.g. iKelloed Il929l : iCohl fc Tohllng [l999l) . 
Although exact, series expansions suffer from a low con- 
vergence rate since these are alternating series; besides, 
the number of integrals increases linearly with the num- 
ber of terms considered up to the truncation order. These 
problems collectively constitute a real practical difhculty 
()Clementl[T97i IStone fc Normanlll992[ ). The proper treat- 
ment of the singularity is the subject of a longstanding chal- 
lenge. Shifting the P-grid and the P'-grid relative to each 



^ This form holds in electrostatics, where p is the density of 
electric charges, and the constant is j^^ (instead of —G). 



other or raising the numerical resolution around the sin- 
gular ity are the most natural techniques (jStemwedel et al.l 
|1990D . but these arc of limited efficiency. When possible, the 
separate treatment of the asymptotic form of the sin gular- 
ity g ives very good results (e.g. lAnsorg et al.l 120031 : iHurel 
l2005f ). although this approach renders the global treat- 
ment somewhat complex. The difficu lty can also be tackled 
by introducing a "softening length" (jHocknev fc EastwoodI 
:198S; Adams et al. 1989). This recipe — widespread in 
disc simulations — must however be seen as nothing but 
a crude approximation that cannot be used for accurate 
modeling under a certain scale. Whatever the prescrip- 
tion for the softening length, which is generally linked 
to the resolutio n or smallest physical length scale (e.g. 
iHure fc PierensI 12009). the inferred force field is globally 
weaker than in the Newtonian case, and the evolution and 
stability of gase ous systems is i nevitably impacted in a non- 
trivial manne r (Romeo 1998: So mmer-Larsen et all 119981 : 
lAdams et anri989:.El-Zant.l998 ). The Poisson equation of 
course provides another othe r way to derive ^ numerically 
(|KelloggJll929l : lDurandlFl953[) . This approach requires accu- 
rate boundary or interior/matching conditions only accessi- 
ble through Eq. ( HI), and complex geometrie s are not always 
easy to manage (|Grandclement et al.ll200ll) . 

In this paper, we present a new means to evaluate Eq.(IT]) 
that avoids the singularity, and, at the same time, properly 
accounts for it. This is achieved by replacing the singular 
kernel by an equivalent and regular, two-term form, one 
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term being the gradient of a new scalar potentiao This is 
the aim of Sections [2] and [3] Kernel equivalence is funda- 
mental to preserving the Newtonian character of the inter- 
action on all scales, and for any separation (in particular 
at long-range). This reformulation is designed to be effi- 
cient within sources, and is not expected to surpass usual 
methods outside sources. From a practical point of view, 
the potential easily becomes accessible as diverging kernels 
have disappeared from the volume integrals. 

Although not specific to a given system of coordi- 
nates, the calculus is developed in cylindrical coordinates 
for which the equivalent kernel takes a nominal form, in 
particular under axial symmetry (this equivalent kernel 
is probably not unique; see Section |4]). It can be applied 
as it is to all rotating gaseous/solid bodies (stars, discs, 
planets, asteroids, etc.) i n either steady s tate or iiot, and 
for various applications jDermottI 119791 : iHachisul Il986bt 
iBaruteau fc Massetl |2008[ ). A basic, six-step algorithm is 
reported in Section [31 together with a numerical experi- 
ment using the most simple quadrature and differentiation 
rules. There is no special assumption about the distribution 
of matter in space (density field and geometry or shape), 
making the method general, and transposable to domains 
of physics other than gravitation. Some interesting perspec- 
tives are listed in the last section. 
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Fig. 1. Typical configuration for the gravitating, celestial 
body, and associated notations. See notelHlfor the definition 
of points P" and Q. 



2. Splitting of the Green kernel 

We consider a volume of space V continuously filled with 
matter, as depicted in Fig.[TJ Using cylindrical coordinates, 
with P'(a, & ^ z) referring to source points, and P(i?, 6*, Z) to 
space points, the above integral for the Newtonian potential 
becomes 



V'(i?,0,z) 



-G 



1 



p{a,9',z)d^T, 



(2) 



^ Some aspects of the t heory of tensor potentials are devel- 
opped by IChandrasekhaii (|1973f ) in the context of rotating, self- 
gravitating fluids. 




Fig. 2. Dimensionlcss Green kernel (left) and cool kernel 
(right) versus a/R for C := and (f) — {0, j, ^, |} labeled 
on the curves. The Green kernel diverges hyperbolically 
when r ^ r' (which corresponds to a — >■ i? and </>—>■ ^ 
here), in contrast to the cool kernel, which remains bounded 
(and even vanishes at the singularity). 



where d?T — adadO' dz is the elementary volume, 

A2 = |r-rf (3) 

= (a + Rf + C^ - 4ai?sin2 (j), 

C = Z-z, (4) 

is the relative altitude, and 

24>^Ti-{e-e'). (5) 

Although A —J- inside V, the potential is generally a 
finite quantity (i.e. integrati on is a regularizing process; e.g. 
lKellogg|[T92alDurandlll95l . We now set 



6 = ^{a + R^ + e > R- (6) 

This quantity is finite everywhere, and non-zero except on 
the polar axis (see below). Assuming i? > 0, we have 

A _ J_A2 



(7) 



1 

A 



1 



4ai?sin^( 
^2 



and so 



1 _ 1 

A" A, 
where we have defined 



1 



+ 4aRsm'q>x-—, (8) 



A(52' 



(9) 



J__ A 

a: = ^- 

The Green kernel is then split into two terms. The term 
■g- is always regulaip ; we call this term the cool kernel in 



^ It is the inverse of a distance, and its value can be interpreted 
geometrically by noting that (see Fig. (TJ 

A 
52 



QP 



/2 



PP 



//2 



d2 



< 00, 



(10) 



where Q{R, 9' + n, Z) is a point of space, diametrically opposite 
to P', and P"(a, 8 + tt,z) is a, point, diametrically opposite to P, 
that belongs to a fictitious source. 
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the following. Figure [3] displays i?/A and -R/A* versus a/R 
around the singularity. We clearly see that the amplitude of 
A^ is bounded, in contrast to A. The second term in Eq.® 
is still singular when A — )■ 0, but its integration over the 
material volume is expected to produce a regular field. The 
idea is to generate this singular kernel from the gradient of 
a regular function k (hereafter called hyperkernel) , and to 
integrate it over the volume V. As the two spaces {R, 9, Z) 
and (a, 9' , z) are independent, the derivative may be drawn 
before the integral. This reasoning can be summarized as 

singular kernel = V regular hyperkernel 



sing, kernel rf'^r = V 



reg. hyperkernel d t, 



where the gradient V is to be taken with respect to one of 
the three variables R, 9, or Z. There are then three possible 
hyperkernels. 

The existence of the hyperkernel is not guaranteed a pri- 
ori, and it is of interest only if it is available in a closed form. 
The investigation indeed shows that the nominal form is ob- 
tained by considering the vertical gradient (i.e. V = dz)- 
This may be due to the special role that the .Z-axis plays 
a in cylindrical coordinates. We therefore do not discuss in 
detail any of the other two options, although these might 
be useful in certain circumstances. 



3. The singular term as the vertical gradient of a 
hyperkernel 

To get the hyperkernel, we consider the integration of the 
singular term in Eq.® with respect to Z. By using the 
intermediate variable t = C/A, we find after some algebra 



dZ 


' f 


S^A 


{a + RYJ 1- 




2^/aR 
m — 



dt 



m? siv? 6 X t'^ ' 



where 



a + R' 
with < TO < 1. We finally get 

dZ 

PA 
' Crn sin < 



K = 4ai? sin 
— m sin (b atanh 



(11) 



(12) 



(13) 



We could obviously add to k, any function of a, i?, 9 and 
0', but this is not necessary here as we take its Z- gradient. 
The Green kernel is then given by the equivalent form 



- = — +dzn. 



(14) 



It is necessary to verify that k is regular. This is straight- 
forward since A > |C| as soon as i? > 0. In other words, if 
m -J> 1 and (/> -^ f , then C,/A -> ±1. 

If we now multiply Eq.® by p(a, 9' , z)d^T — which does 
not depend on Z — and integrate over the material volume 
V, we find that 



V ^ 






dz 



pnd^T, (15) 



where the partial derivative now operates on the integral. 
Up to a factor — G, this expression is precisely the poten- 
tial defined by Eq.([T]), and it is both exact and general. 
It depends on neither the body's shape nor on the distri- 
bution of its mass density. It applies not only to volume 
distributions, but also to surface distributions (see Sect. H]) 
and linear distributions. The first integral in Eg. dTS]) is then 
the cool potential associated with the cool kernel, and the 
second term is the vertical gradient of a hyperpotential. 

On the polar axis (i.e. R = 0) A = 6, and so Eq.® 
does not help us to treat the singularity when C = and 
a = 0. In this case, we have 



K — -—dZ 
J A 

= asinh — , 
a 



(16) 



provided that a > 0. The potential can then be written 



^iO,Z) 



-Gdz 



pnd^ 



(17) 



where here there is no cool kernel. Besides, we see that 
liuia-^o aK = 0. 

There is no continuity between the two different expres- 
sions for the cool kernel, the one valid at i? = and the 
other valid at _R — > (this is also true for the hyperkernel) . 
This is no problem as long as we do not have to consider the 
radial gradient of n (see below the numerical experiment). 

Finally, we note that, as we work with an equiva- 
lent form of the Green kernel, the potential found from 
Ea.([T5|) automatically has the right asymptotic property, 
and varies like M/r sufficiently far away from the body. At 
large relative distance (i.e. R ^ a and Z ^ z), we have 
S 



= \/i?2 + Z^. At the lowest order, one finds that 



1 j3 ^I 



(18) 



which means that the long-range behavior is exclusively 
contained in the cool kernel. 



4. The case of axially-symmetric bodies 

Axially-symmetric bodies constitute an important class 
of ast rophysical objects (jChandrasekhad Il973t iHachisul 
Il986al ). Interestingly enough, in problems where do'P = 0, 
we can rewrite the above expressions in a more compact 
form in terms of elliptic integrals. The first integral in the 
right-hand-side of Eq. p3|) becomeQ 




—kE{k)dadz, (19) 
R 



where 



E(fc) 



7r/2 



V 1 — fc^ sin X dx, 



(20) 



is the complete elliptic integral of the second kind, 
k — 2^/aR/5 is the modulus (with k e [0, 1]), and the 



'* A factor of two is due to dO' /d(f>, and another factor of two 
contained in the modulus k comes from symmetry considera- 
tion (i.e. matter located at 9' G [9,6 + n] provides the same 
contribution as matter located at 9' £ \9,9 ~ n]). 
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double integration runs over the meridional cross section S 
of the body. The integration over the polar angle 9' of the 
hyperkernel gives: 



sin (h atanh 



^m sin ( 



cos 6 atanh 



^m sin ( 



+ 



mS 



F{(l),k)-m' n{(j),m,k) 



(21) 



where 



i^(</>,fc) 



dx 



'0 v 1 — k'^ sin^ X 
is the incomplete elliptic integral of the first kind, and 

/"^ dx 

n((/), m, k) 



(l — jTi^ sin x) vl — k^ 



■ 2 
sm X 



'221 



(23) 



is the incomplete elliptic integral of the third kind (to is the 
parameter and m' = \J\ — m?). Over the whole circle (i.e. 
(j) G [0, f ]), this yields the axially symmetric potentiao 



ipiR,Z) 



2G 11 pJ -kF,(k)dadz 




(25) 



P\l — CH(to., k)dadz 
5 V it 



where H is defined for convenience by 

H(?7i, fc) = fc K(A:) — ?7i n(?7i, fc) 



(26) 



with n(TO, k) = n(|, TO, k) and K(fc) = F{^,k). The pres- 
ence of the H-functi on is actually expe cted here since dzA 
and dzA are linked f|Trova et alil2012t) . 
On the polar axis, we get 



^"(0, Z) = —2'iTGdz I / pnadadz, 



(27) 



where k is, in this case, given by Ea. (jl6p . We note that, in 
this axially symmetric case, the potential could be deter- 
mined through Eg. p^ (i.e. by using the cool kernel only, 
and no hyperkernel), but the integrand still contains a hy- 
perbolic divergence as a — >■ and C —?' which is not easy 
to manage. This is why it seems much better to consider 
Ea. (j27p . as the logarithmic divergence of the hyperkernel 
(i.e. the asinh term) is cancelled out by the elementary 
volume when a is close to (i.e. lima_i.o an ~ 0). 



1. Discretize the material source V on a 3D-grid, in the 
form of quadruplets {P'(aij^fe, ^ij,fe, ^jj.fc), Pij,fc}, 

2. Select a point P{R, 9, Z) in space, 

3. For each point P'ij.fe of the source, compute the cool 
kernel, and the hyperkernel k from Eqs.® and (fT3)) or 
Eq.llll), 

4. Perform the volume integrals in Ea. p^ or (|T7| . 

5. Determine the vertical gradient of the hyperpotential, 

6. Add this derivative to the cool potential, and multiply 
by-G, 

and reiterate steps 2 to 6 to generate a potential map. The 
components of the associated gravitational force are de- 
duced as usual from the three gradients of V'- 

We have considered a homogeneous, axially symmet- 
ric torus with a square meridional cross section S with 



(a, z) e 



[1 31 



2' 2J 



The mass density {p — 1 inside 



S and outside) is defined on a regular mesh, with N 
nodes in each direction. The computational grid also con- 
sists of a square box (i?, Z) € [0, 2] x [—1, 1] with L nodes 
per direction and regular spacing, therefore encompassing 
the body. The double integrals in Eas. ((25)l and ((27l) are 
computed at each node of the computational grid through 
the two-point trapezoidal rule. The partial derivative of the 
hyperpotential is estimated through second-order finite dif- 
ferences. These basic schemes are very easy to implement 
and the above six-step procedure contains no pitfall. We 
take iV = L = 31 here. Figures [3^-d display, respectively, 
the integral of the cool kernel, the integral of the hyperker- 
nel K, its vertical gradient, and the total potential, obtained 
by adding the first and third maps. The boundary of the 
toroidal body is superimposed on these maps. We note that 
the vertical gradient of the hyperpotential makes the geo- 
metrical cross section rise above the background, and filters 
the curvature effects (which are enhanced by the integral 
of the cool kernel). 

The computing time is typical of integral methods. 
Under the conditions of the present example, the integra- 
tion of the cool kernel and hyperkernel (step 4 above) re- 
quires 2N'^ elementary operations per point of space. To get 
the potential at the nodes oi a L x L square grid, we then 
find 2N'^L'^ (the time needed to determine the vertical gra- 
dients is negligible in comparison). This is therefore much 
smaller than that is usually obtained from an expanded 
Green function by a factor equal to the number of terms 
up to the truncation order (which can be as high as a few 
hundred; see e.g. iHachisuil986a : iStone &: Norma n 1992.) . 



5. An example of implementation 

We now present a first numerical experiment to briefly de- 
scribe the main steps of the method, and demonstrate its 
simple implementation. For a full 3D body, the main stepqj 
of any numerical estimate can be summarized as 



^ If we perform the Z-derivati ve and rearran ge terms, we re- 
cover the well-known expression (|Durandll 19531 ). namely 



^{R,Z) = 



-kK.[k)dadz, 



(24) 



whose kernel is logarithmically singular. 

^ Step 1 should be executed once and for all, except for sys- 
tems that evolve with time. 



6. Checking the accuracy 

The numerical accuracy is sensitive to various ingredients, 
such as the quadrature and differentiation schemes. To a 
lesser extent, it also depends on the mass density distribu- 
tion p and equation of the boundary 9V, which may gen- 
erate additional difficulties in the calculations (interpola- 
tion of data points, infinite derivatives at edges, etc.). We 
present a second numerical test illustrating the accuracy 
of the method by considering an exact potential/ density 
pair. Not many configur ations correspond to finite m ass and 
finite size systems fe.g. iBinnev &: Tremainelll987t ). When 
matter is gathered in a plane (i.e. a flat disc), Ea.([T5| is di- 
rectly transposable by setting p{a, 9', z) = S(a, 9')5{z) and 
integrating over z. Under axial symmetry, we respectively 
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Fig. 3. Main steps in the computation of the potential from Eas. (l25]) and (|27)) in a typical case: (a) the integral of the 
cool kernel, (b) the integral of the hyperkernel, (c) its vertical gradient, and (d) the potential as the sum of maps (a) 
and (c). Here, the body is an axially symmetric torus with square cross section (boundary indicated with a white line). 
As clearly visible in graphs (a)-(c), the treatment differs slightly on the polar axis (first column of pixels). See the text 
for the numerical setup. 



ged from Eqs.dlSl) and ([27] 



where a-m and Cout denote the discs inner and outer edges, 
and 



iP{R,Z) 



2G 



2Gdz 



— koE{ko)da 
it 



Si / — ZH(m, fco)da 
R 



(29) 



P - 

Air, 



4aR 



{a + i?)2 + Z2 ■ 



(31) 



for R> 0, and 



^iO,Z) 



-2TTGdz 



Yi asinh 



ada. 



"^ Thes e two formulae can be compared with th e classical ex- 
pression (|Durandlll95^ : lBinnev fc Tremainelll987n : 



ip{R,Z) = -2G 



T., —ko'K{ko)da, 
H 



which is singular inside the source. 



By setting ain = and E ex (1 — a^/oo^J^/^ for g g 
[0, «ou t] , one gets one of the three cases analyzed by ISchulzj 
(|2009[ ). For such a distribution, the associated potential, 

(30) ipe, is known exactly in a closed-form for any point of space. 
Figure m gives tpe as well as the error index e = log^g |1 — 
ip/ipel where V' is determined from Eas. ip^ and ([30)1 in 
the same conditions as above (we set Cout = !)■ We see 
that the potential outside and especially inside the disc is 
well-reproduced. The relative error, on the order of 10"'^, 

(28) agrees with the second-order of the schemes at the actual 
mesh size of ""j^"""' = ■^. The accuracy can be tuned by 
changing the quadrature and differentation schemes. 
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Fig. 4. Same legend as for Fig. |3]but for the flat, axially symmetric disc with surface density E ex (1 



,2N3/2 



inner edge 



Oin = and outer edgeoout = 1- The disc is indicated with a white line. The exact potential ipe (left) is derived from the 
formula of Schuld (|2009r ). The associated error index (right) is determined once i/j is computed from Eas. ((29t and ((30)) . 
The mesh size is -^ corresponding to A^ = 31 radial points. 



7. Concluding remarks 

We have reformulated the Green kernel appearing in po- 
tential problems to circumvent the singularity and, at the 
same time, properly account for it. As a consequence, the 
gravitational potential of any celestial body, regardless its 
shape and matter density distribution, becomes directly ac- 
cessible through two "classical" volume integrals, followed 
by a partial derivative. The method is applicable to three- 
dimensional, fully inhomogeneous systems, as well as to sur- 
face and line distributions. It is especially efficient inside 
distributions where most approaches exhibit a real prac- 
tical complexity, converge very slowly, or produce spuri- 
ous errors.. The presence of regular kernels ensures that 
the method is stable and easy to implement. This should 
encourage modellers to abandon various integration tech- 
niques that do not "faithfully" reproduce the Newtonian 
character of the potential and forces. In the context of discs 
for instance, this method appears to be a real alternative 
to softened gravity, which remains a free parameter, non- 
Newtonian theory. As stressed, it is probably possible to 
determine other cool kernel/hyperkernel pairs (for instance, 
by considering the hyperkernel as a radial/angular gradi- 
ent), but the one presented in the body of this paper seems 
the simplest one. It is in particular well-suited to axially 
symmetric configurations. 

This study needs to be continued in several respects, in- 
cluding the analysis of the mathematical properties of the 
cool kernel and hyperkernel and their physical meanings, as 
well as the derivation of the equivalent kernel in other sys- 
tems of coordinates (e.g. cartesian and spherical coordinate 
systems) . In addition, it would be interesting to expand the 
two kernels in Ea. (jl5|) in series, and compare their proper- 
ties with the expansion of the Green function in Legendre 
polynomials, inside as well as outside the body. The cool 
kernel/hyperkernel pair is also interesting as a new start- 
ing point to generating various kinds of approximations. 
Apart from the astrophysical context where there are so 



many applications about gravitation, this technique is also 
transposable to other kinds of problems involving improper 
integrals. This can be, for instance in electromagnetism, 
the determination of the potential vector A and associated 
magnetic field induce d by current densities (iJac kson 1992 
'Cohl k Tohlinel ll999D . These points will be touched on in 
forthcoming papers. 
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